import geopandas as gpd
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import networkx as nx
import shapely
import copy, os, sys
from scipy.interpolate import interp1d
from scipy.signal import savgol_filter
from scipy.misc import derivative
# from typing_extensions import final#求导用
sys.path.append(os.path.abspath("../"))
from merge.utils import move_extend
from utils.basic_function import *
from settings import *

plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False


class MethodJoinCut:
    """
    融合处理工具类，对数据集对象内的上报路网与遥感路网进行融合
    """
    def __init__(self,
                 dataprovider,
                 THRESHOLD_BUFFER_COVER=0.6,
                 THRESHOLD_TORLERANCE=30
                 ):
        """初始化变化检测融合处理工具对象

        Args:
            dataprovider (DataSetProvider)): 待处理数据集对象
            THRESHOLD_BUFFER_COVER (float, optional): 判断裁剪后的路段与上报路网相似程度阈值的因子. Defaults to 0.6.
            THRESHOLD_TORLERANCE (int, optional): 相邻线段捕捉距离阈值，仅在使用shapely方法与fix方法时有效. Defaults to 30.
        """
        self.dataprovider = dataprovider
        self.match_df = dataprovider.match_result
        self.rp_df = dataprovider.rp_df
        self.THRESHOLD_BUFFER_COVER = THRESHOLD_BUFFER_COVER
        self.THRESHOLD_TORLERANCE = THRESHOLD_TORLERANCE
        self.SRID = dataprovider.SRID
        pass

    def my_extend(self, geom_self, geom_base, tolerance=30):
        """
        改进的捕捉工具，对应fix方法
        把geom_self的端点捕捉到geom_base上。
        · 使用geom_self离geom_base最近的端点
        · 若geom_self的端点离geom_base的端点的距离符合要求，优先捕捉到geom_base的端点上
        · 若不符合要求，捕捉到离geom_self端点最近的点上。

        Args:
            geom_self ([type]): 当前捕捉对象
            geom_base ([type]): 目标捕捉对象
            tolerance (int, optional):捕捉距离阈值. Defaults to 30.

        Raises:
            ExtendException: 捕捉错误

        Returns:
            shapely.geometry: 返回捕捉后的地理对象
        """

        st_self = geom_self.coords[0]
        ed_self = geom_self.coords[-1]
        st_base = geom_base.coords[0]
        ed_base = geom_base.coords[-1]
        st_dist = geom_base.distance(Point(st_self))
        ed_dist = geom_base.distance(Point(ed_self))
        if st_dist <= ed_dist and st_dist < tolerance:
            # 如果端点也符合条件，捕捉到端点上
            if cal_dist(st_self, st_base) < tolerance:
                pt_proj = st_base
            elif cal_dist(st_self, ed_base) < tolerance:
                pt_proj = ed_base
            else:
                pt_proj = shapely.ops.nearest_points(Point(st_self),
                                                     geom_base)[1].coords[0]
            list_new = [pt_proj] + list(geom_self.coords)[1:]
            assert cal_dist(st_self, pt_proj) < tolerance
        elif ed_dist < st_dist and ed_dist < tolerance:
            # 如果端点也符合条件，捕捉到端点上
            if cal_dist(ed_self, st_base) < tolerance:
                pt_proj = st_base
            elif cal_dist(ed_self, ed_base) < tolerance:
                pt_proj = ed_base
            else:
                pt_proj = shapely.ops.earest_points(Point(ed_self),
                                                    geom_base)[1].coords[0]
            list_new = list(geom_self.coords)[:-1] + [pt_proj]
            assert cal_dist(ed_self, pt_proj) < tolerance
        else:
            raise ExtendException(min(st_dist, ed_dist), tolerance)
        geom_self_new = LineString(list_new)
        return geom_self_new

    def snap_fix(self, geo1, geo2, tolerance=20):
        """改进的捕捉策略，已弃用

        Args:
            geo1 ([type]): [description]
            geo2 ([type]): [description]
            tolerance (int, optional): [description]. Defaults to 20.

        Returns:
            [type]: [description]
        """
        g11 = np.array(geo1.coords[0])
        g12 = np.array(geo1.coords[-1])
        g21 = np.array(geo2.coords[0])
        g22 = np.array(geo2.coords[-1])
        g1 = [g11, g12]
        g2 = [g21, g22]
        ax = np.array([
            np.sum((np.array(g21) - np.array([g11, g12]))**2, axis=1),
            np.sum((np.array(g22) - np.array([g11, g12]))**2, axis=1)
        ])
        # (np.array(g21) - np.array([g11,g12]))**2
        index = np.argmin(ax)
        index1 = int(index / 2)
        index2 = int(index % 2)
        ## 这样可以确定最近的两个点，然后连接起来 再merge

        join_geo = shapely.geometry.LineString([g2[index1], g1[index2]])
        if (join_geo.length > tolerance):
            return geo2
        else:
            return shapely.ops.linemerge([join_geo, geo1])

    ## 剪断路网
    def cut(self, line_series, distance: float, NID=0) -> pd.DataFrame:
        """
        在距离线段起点某一指定距离处打断路段

        Args:
            line_series (geopandas.GeoSeries): 源路段对象
            distance (float): 距离路段起点长度
            NID (int, optional): 打断点编号. Defaults to 0.

        Returns:
            pd.DataFrame: 打断后的路段集合
        """
        line = line_series.geom

        ## * 打断点在起点前
        if distance <= 0.0:
            line_cut_series1 = pd.Series({
                "geom":
                shapely.geometry.LineString(line),
                'start':
                NID,
                'end':
                line_series.end
            })
            # ldf = pd.DataFrame([shapely.geometry.LineString(line),NID,line_series.end]).T  #-1 end -2 start
            ldf = pd.DataFrame([line_cut_series1])
            ldf.columns = ['geom', 'start', 'end']
            return ldf

        ## * 打断点在起点后
        elif distance >= line.length:
            line_cut_series1 = pd.Series({
                "geom":
                shapely.geometry.LineString(line),
                'start':
                line_series.start,
                'end':
                NID
            })
            # ldf = pd.DataFrame([shapely.geometry.LineString(line),line_series.start,NID]).T  #-1 end -2 start
            ldf = pd.DataFrame([line_cut_series1])
            ldf.columns = ['geom', 'start', 'end']
            return ldf

        ## * 打断点在路段中
        coords = list(line.coords)
        ### 循环遍历地理图形内部结点
        for i, p in enumerate(coords):
            project_distance = line.project(shapely.geometry.Point(p))
            ## 断点刚好位于内部结点上
            if project_distance == distance:
                line_cut_series1 = pd.Series({
                    "geom":
                    shapely.geometry.LineString(coords[i:]),
                    'start':
                    NID,
                    'end':
                    line_series.end
                })
                line_cut_series2 = pd.Series({
                    "geom":
                    shapely.geometry.LineString(coords[:i + 1]),
                    'start':
                    line_series.start,
                    'end':
                    NID
                })
                ldf = pd.DataFrame([line_cut_series1, line_cut_series2])
                # ldf = pd.DataFrame([
                #     [shapely.geometry.LineString(coords[i:]),shapely.geometry.LineString(coords[:i+1])],
                #     [NID,line_series.end],
                #     [line_series.start,NID]
                #     ]).T  #-1 end -2 start
                # ldf.columns = ['lines','start','end']
                return ldf
            ## 打断点位于两节点中间
            if project_distance > distance:
                cp = line.interpolate(distance)
                line_cut_series1 = pd.Series({
                    "geom":
                    shapely.geometry.LineString([(cp.x, cp.y)] + coords[i:]),
                    'start':
                    NID,
                    'end':
                    line_series.end
                })
                line_cut_series2 = pd.Series({
                    "geom":
                    shapely.geometry.LineString(coords[:i] + [(cp.x, cp.y)]),
                    'start':
                    line_series.start,
                    'end':
                    NID
                })
                ldf = pd.DataFrame([line_cut_series1, line_cut_series2])

                # ldf = pd.DataFrame([
                #     [shapely.geometry.LineString([(cp.x, cp.y)] + coords[i:]), shapely.geometry.LineString(coords[:i] + [(cp.x, cp.y)])],
                #     [NID,line_series.end],
                #     [line_series.start,NID]
                #     ]).T  #-1 end -2 start
                ldf.columns = ['geom', 'start', 'end']
                return ldf

                # return [
                #     shapely.geometry.LineString([(cp.x, cp.y)] + coords[i:]),
                #     shapely.geometry.LineString(coords[:i] + [(cp.x, cp.y)])
                #     ]
        pass

    def cut_by_multidistance(self, line, distances, NIDs=None):
        """用一组描述打断点到路段起点的距离集合打断指定路段

        Args:
            line ([type]): 源路段图形对象
            distances (list): 打断点到路段起点的距离集合 MUST BE a LIST which is SORTED by distace DESC
            NIDs ([type], optional): 打断点编号集. Defaults to None.

        Returns:
            geopandas.GeoDataFrame: 打断路段集合
        """
        if (NIDs is None):
            NIDs = range(len(distances))
        if (len(distances) == 0):
            return [line]
        line_series = pd.Series({"geom": line, 'start': -2, 'end': -1})
        ## ** DOing
        # distances.sort(reverse=True)
        lineDf = self.cut(line_series, distances[0], NIDs[0])
        for i in range(1, len(distances)):
            distance = distances[i]
            NID = NIDs[i]
            cuttingline = lineDf.iloc[-1]
            lineDf.drop([len(lineDf) - 1], inplace=True)
            lineDf = lineDf.append(self.cut(cuttingline, distance, NID))
            lineDf.reset_index(drop=True, inplace=True)
        return lineDf

    def cut_by_multipoints(self, pois_df, geom):
        """
        在多点处打断目标路段

        Args:
            pois_df (geopandas.GeoDataFrame): 打断点集合
            geom (shapely.geometry]): 源路段图形对象

        Returns:
            geopandas.GeoDataFrame: 打断路段集合
        """
        distances = []
        for poi_i in range(len(pois_df)):
            poi = pois_df.iloc[poi_i].pois
            distances.append(geom.project(shapely.geometry.Point(poi)))
        pois_df["distance"] = distances
        pois_df.sort_values(by='distance', inplace=True, ascending=False)
        cutlines = self.cut_by_multidistance(geom, pois_df.distance.values,
                                             pois_df.NID.values)
        return cutlines

    ## 获取距离分布函数
    def getHist(self, rp_geom, rs_geom, step=5, density_method='manually'):
        """获取上报路段-遥感路段距离变化函数

        Args:
            rp_geom (shapely.geometry): 上报路段图形对象
            rs_geom (shapely.geometry): 遥感路段图形对象
            step (int, optional): 程序自动增密步长. Defaults to 5.
            density_method (str, optional): 增密方法，默认程序内手动实现增密，否则不增密. Defaults to 'manually'.

        Returns:
            list: [遥感路网上计算点到上报路网的距离集，两距离变化值，插入点在线性坐标系中的坐标，计算点集]
        """
        # tmp_match_rodas_df = match_df_todo[match_df_todo.MATCH_ROAD == ldbm]
        # tmp_rp_rodas_df = rp_c_df[rp_c_df.ldbm == ldbm]
        geom = rp_geom
        # geom = tmp_rp_rodas_df.geom.values[0]
        # pois = list(tmp_match_rodas_df.geometry.values[0].coords)
        # for i in range(1,len(tmp_match_rodas_df.geometry.values)):
        #     pois.extend(list(tmp_match_rodas_df.geometry.values[i].coords))
        # * 一种思路是使用QGIS的加密功能呢先进行加密后计算 缺点是，点分布不均匀，但是重要的特征点都能保留 -------好像也没什么影响？ 20210715
        # ? 又想了一下，好像是有影响的，太密的地方会出问题
        # * 手动加密一下
        if (density_method == 'manually'):
            pois = [
                rs_geom.interpolate(inter_dist)
                for inter_dist in range(0,
                                        int(rs_geom.length) + 2, step)
            ]
        else:
            pois = list(rs_geom.coords)
        dists = []
        cut_i = -1
        interpolate_dists = []
        for i, poi in enumerate(pois):
            src_poi = shapely.geometry.Point(poi)
            dist = geom.project(src_poi)
            interpolate_dists.append(dist)
            proj_poi = geom.interpolate(dist)
            dp = proj_poi.distance(src_poi)
            dists.append(dp)
        dists = np.array(dists)
        ddist = dists[1:] - dists[0:-1]
        return dists, ddist, interpolate_dists, pois

    ## 路段合并函数
    ## 先合并了  但是把节点提取出来了
    def joinSegments(self, geoms):
        """合并多个路段

        Args:
            geoms (list)): 合并路段集

        Returns:
            list: [合并点，合并后路段图形]
        """
        joinLine = shapely.ops.linemerge(list(geoms))
        pois = []
        g = nx.DiGraph()
        for i in range(len(geoms)):
            # sp = shapely.geometry.Point(*match_df.iloc[i].geometry.coords[0])
            # ep = shapely.geometry.Point(*match_df.iloc[i].geometry.coords[-1])
            # pois.append(sp)
            # pois.append(ep)
            g.add_edge(geoms[i].coords[0], geoms[i].coords[-1])
        for poi in list(g.degree):
            if (poi[-1] > 1):
                pois.append(shapely.geometry.Point(poi[0]))
        return pois, joinLine

    ## 线性内插插入点距离
    def interpolateFromX(self, x: np.array, interpolate_dists: list):
        return ((np.array(interpolate_dists)[np.ceil(x).astype(int)] - np.array(interpolate_dists)[np.floor(x).astype(int)]) * \
        (x - np.floor(x) ) + np.array(interpolate_dists)[np.floor(x).astype(int)]).reshape(-1).tolist()

    def get_cutting_position(
        self,
        ldbm,
        dists,
        interpolate_dists,
        distance_threshold=MethodjoinCutThresholdSettings.
        distance_threshold,
        filter_window_length=MethodjoinCutThresholdSettings.
        filter_window_length,
        x_length=MethodjoinCutThresholdSettings.x_length,
        derivative_threshold=MethodjoinCutThresholdSettings.
        derivative_threshold,
        filter_window__min_length=MethodjoinCutThresholdSettings.
        filter_window__min_length,
        savgol_filter_polyorder=MethodjoinCutThresholdSettings.
        savgol_filter_polyorder):  ## 还有好多threshold 、 plot 明天设置吧 好困啊
        """获取打断点位置

        Args:
            ldbm ([type]): 待打断路段的路段编码
            dists ([type]): 距离变化值
            interpolate_dists ([type]): 计算点线性坐标值
            distance_threshold ([type], optional): 变化阈值. Defaults to Method_joinCut_Threshold_Settings.distance_threshold 为30.
            filter_window_length ([type], optional): 滤波窗口大小. Defaults to Method_joinCut_Threshold_Settings.filter_window_length 为27.
            x_length ([type], optional): 求导细粒度. Defaults to Method_joinCut_Threshold_Settings.x_length.
            derivative_threshold ([type], optional): 导数阈值，小于该阈值认为是极值点. Defaults to Method_joinCut_Threshold_Settings.derivative_threshold 为0.5.
            filter_window__min_length ([type], optional): 最小滤波窗口大小. Defaults to Method_joinCut_Threshold_Settings.filter_window__min_length 为5.
            savgol_filter_polyorder ([type], optional): [description]. Defaults to Method_joinCut_Threshold_Settings.savgol_filter_polyorder.

        Returns:
            list: 打断点线性坐标，打断点原始数据
        """

        if (len(dists) < filter_window_length):
            filter_window_length = filter_window__min_length
        filter_dists = savgol_filter(dists, filter_window_length,
                                     savgol_filter_polyorder)
        f_filter = interp1d(range(len(dists)),
                            filter_dists,
                            kind="cubic",
                            fill_value="extrapolate")
        x = np.linspace(
            1,
            len(dists) - 2,
            x_length)  ## ! 会遇到超限的问题，需要可以外推 fill_value="extrapolate"
        y_filter = f_filter(x)

        ## 求导
        a = np.empty((len(x)))
        for i in range(len(a)):
            a[i] = derivative(f_filter, x[i],
                              x[1] - x[0])  #求导，第二个参数是求导位置，第三个参数是dx
        ## 极值点
        index = np.argwhere(abs(a) < derivative_threshold).reshape(-1)
        ## 加入两端
        index2 = np.hstack([np.array([0]), index, np.array([len(x) - 1])])

        ## 筛选无关点
        delta = y_filter[index2[1:]] - y_filter[index2[0:-1]]
        index_delta_positive = np.argwhere(delta > distance_threshold)
        index_delta_negative = np.argwhere(delta < -distance_threshold)

        ## 找出极值点所在的index
        x_positive = x[index2[0:-1][index_delta_positive]].reshape(-1)
        x_negative = x[index2[1:][index_delta_negative]].reshape(-1)

        x_map = []
        print("length x_po", len(x_positive))
        print("length x_ne", len(x_negative))

        for _x in x_positive:
            x_map.append([_x, 1])
            pass
        for _x in x_negative:
            x_map.append([_x, 2])
            pass
        print("length x_map", x_map)
        # x_map = sorted(x_map.items(), key=lambda d: d[0])
        x_map.sort(key=lambda d: d[0])
        print("length x_map", x_map)
        ## 在内插距离表中内插极值点在上报路段中的线性系统坐标
        cutDists = self.interpolateFromX(x_positive,interpolate_dists) + \
                self.interpolateFromX(x_negative,interpolate_dists)

        # *绘图
        if (image_save_option):
            # # ! 临时施工
            fig, axes = plt.subplots(figsize=(10, 5), nrows=1, ncols=1)
            axes.plot(x, y_filter, label="平滑后")
            axes.plot(range(len(dists)), dists, alpha=0.5, label="平滑前")
            axes.legend(fontsize=20)
            axes.set_title("距离分布函数", fontsize=20)
            # plt.show()

            fig, axes = plt.subplots(figsize=(6, 12), nrows=3, ncols=1)
            axes[0].scatter(x[index], a[index], label="极值点候选点")
            axes[0].plot([0, len(dists)], [0, 0])
            axes[0].plot(x, a, label="导数")
            axes[0].set_title("导数")
            axes[1].set_title("候选点位置")
            axes[1].plot(x, y_filter, label="距离函数")
            axes[1].scatter(x[index], y_filter[index])
            axes[2].plot(x, y_filter)
            axes[2].scatter(x_positive,
                            y_filter[index2[0:-1][index_delta_positive]],
                            label="正裁剪点")
            axes[2].scatter(x_negative,
                            y_filter[index2[1:][index_delta_negative]],
                            label="负裁剪点")
            axes[2].set_title("裁剪点位置")
            axes[2].legend()
            # plt.show()

            # # *----------------------------------------------------------
            fig, axes = plt.subplots(figsize=(10, 5), nrows=1, ncols=4)
            ax1 = axes[0]
            ax2 = axes[1]
            ax3 = axes[2]
            ax4 = axes[3]
            ax1.plot(x, y_filter)
            ax1.plot(range(len(dists)), dists, alpha=0.5)
            ax1.set_title("origin & sg_filter")

            ax2.scatter(x[index], a[index])
            ax2.plot([0, len(dists)], [0, 0])
            ax2.plot(x, a)
            ax2.set_title("dx")

            ax3.plot(x, y_filter, label="filter")
            ax3.scatter(x[index], y_filter[index])
            ax3.set_title("position")

            ax4.plot(x, y_filter)
            ax4.scatter(x_positive,
                        y_filter[index2[0:-1][index_delta_positive]])
            ax4.scatter(x_negative, y_filter[index2[1:][index_delta_negative]])
            ax4.set_title("positive & negetive Cutting position")
            plt.suptitle(ldbm)
            # self.viz.matplot(plt)
            plt.savefig(f"{image_path}/{ldbm}_CuttingPois.png")

        return cutDists, x_map

    def getPoisFromDists(self, dists: list, geom):
        pois = []
        for dist in dists:
            poi = geom.interpolate(dist)
            pois.append(poi)
        return pois

    def joinMatchCutWorkflow_single(self,
                                    ldbm,
                                    geom_rs,
                                    geom_rp,
                                    joinMethod='shapely'):
        """合并单条上报路网及其相关遥感路网工作流程
        思路：
        1. 寻找切割点
        2. 切割线段
        3. 计算相似程度与合并分配
        4. 合并路网
        5. 连接与捕捉

        Args:
            ldbm (str): [description]
            match_df (geopandas.GeoDataFrame): Roads from remote sensing
            rp_df (geopandas.GeoDataFrame): Roads from reporting data
            THRESHOLD_TORLERANCE (int, optional): snap threshold. Defaults to 50.
            joinMethod (str, optional): method to snap. 
                                        1."shapely"-->the method from shapely,but it seems that it is not always right for some strange Roads in some situation .
                                        2."fix" --> the method which is simply snap the vertices of the roads from rs to road from rp,but it performs not so well at the snap poi. 
                                        Defaults to 'shapely'.
        """
        print(f"****************************正在处理{ldbm}************************")
        if (len(geom_rs) > 1):
            print("路段数：", len(geom_rs))
        segment_pois, geom_rs = self.joinSegments(geom_rs)
        linePois = list(geom_rs.coords)

        ## * 1. 寻找裁剪点
        dists, ddist, interpolate_dists, pois = self.getHist(geom_rp, geom_rs)
        cutting_dists, x_map = self.get_cutting_position(
            ldbm, dists, interpolate_dists, filter_window_length=27)
        pois = self.getPoisFromDists(cutting_dists, geom_rp)
        poi_df = pd.DataFrame([pois, ["cut"] * len(pois)]).T
        poi_df.columns = ['pois', 'kind']
        poi_df_temp = pd.DataFrame([segment_pois,
                                    ["seg"] * len(segment_pois)]).T
        poi_df_temp.columns = ['pois', 'kind']
        poi_df = poi_df.append(poi_df_temp)

        # pois.extend([shapely.geometry.Point(geom_rp.coords[0]),shapely.geometry.Point(geom_rp.coords[-1])]) ## 首尾点
        ## * 2.检查首尾点是否需要加入裁剪
        head_tail_pois = [
            shapely.geometry.Point(*linePois[0]),
            shapely.geometry.Point(*linePois[-1])
        ]
        ### 如果是空 意味着中间没有裁剪 所以是匹配上的，应放入前后
        if (len(cutting_dists) == 0):
            poi_df_temp = pd.DataFrame(
                [head_tail_pois, ["headtail"] * len(head_tail_pois)]).T
            poi_df_temp.columns = ['pois', 'kind']
            poi_df = poi_df.append(poi_df_temp)
            pois.extend([
                shapely.geometry.Point(*linePois[0]),
                shapely.geometry.Point(*linePois[-1])
            ])  ## 首尾点
        ### 如果裁剪点不空，中间裁剪，两端还需要判断是否需要
        else:
            # 2表示negative 1表示positive
            if (x_map[0][1] == 2 and dists[0] > 30):  ## TODO这里还有个条件
                head_tail_pois.pop(0)
            if (x_map[-1][1] == 1 and dists[-1] > 30):
                head_tail_pois.pop(-1)
            poi_df_temp = pd.DataFrame(
                [head_tail_pois, ["headtail"] * len(head_tail_pois)]).T
            poi_df_temp.columns = ['pois', 'kind']
            poi_df = poi_df.append(poi_df_temp)
            pois.extend(head_tail_pois)  ## 首尾点
            pass
        pass

        ### TODO 整理poi_df 目的是为了删除一些离起点终点很近的打断点，阈值为20 这里有个阈值 20210818
        # poi_df中的所有点都是要参与打断路段的，以下步骤是整理这些点，排除掉离得比较近的点，减少错误
        poi_df["NID"] = range(len(poi_df))
        poi_df.reset_index(drop=True, inplace=True)

        poi_seg_df = poi_df[poi_df.kind != 'cut']
        poi_i = 0 # 作为指针指向待处理的数据行
        while (poi_i < len(poi_df)):
            poi = poi_df.iloc[poi_i]
            if (poi.kind != 'cut'):
                ## 只处理裁剪点
                poi_i += 1
                continue
            deleteTag = False
            poi_cut = poi_df.iloc[poi_i].pois
            # 检查当前裁剪点与其他非裁剪点之间是否过近
            for poi_seg_i in range(len(poi_seg_df)):
                poi_seg = poi_seg_df.iloc[poi_seg_i].pois
                if (poi_cut.distance(poi_seg)) < 20:
                    poi_df.drop(index=[poi_i], inplace=True)
                    poi_df.reset_index(drop=True, inplace=True)
                    deleteTag = True
                    break
            if (not deleteTag):
                poi_i += 1
            pass
        pass
        ## * 3.切割线段
        cutlines_rp = self.cut_by_multipoints(poi_df, geom_rp)
        cutlines_rs = self.cut_by_multipoints(poi_df, geom_rs)

        df = pd.DataFrame(cutlines_rp)
        df["NID"] = range(len(df))
        df["NID_rs"] = [-1] *len(df)
        df = gpd.GeoDataFrame(df, geometry="geom")

        df_rs = pd.DataFrame(cutlines_rs)
        df_rs["NID"] = range(len(df_rs))
        df_rs = gpd.GeoDataFrame(df_rs, geometry="geom")

        ## * 4. 计算切割后图形的相似程度
        hausdorff_distance = np.ones([len(cutlines_rs),
                                      len(cutlines_rp)]) * 1e6
        bufferPer = np.ones([len(cutlines_rs), len(cutlines_rp)]) * 1e6
        for i in range(len(cutlines_rs)):
            geo = df_rs.iloc[i].geom
            for j in range(len(cutlines_rp)):
                cut_rp_geo = df.iloc[j].geom
                hdist = geo.hausdorff_distance(cut_rp_geo)
                hausdorff_distance[i, j] = hdist
                ## TODO 还有个阈值在这
                bufferPer[i, j] = geo.buffer(
                    MethodjoinCutThresholdSettings.check_similarity_buffer
                ).intersection(
                    cut_rp_geo.buffer(
                        MethodjoinCutThresholdSettings.
                        check_similarity_buffer)).area / geo.buffer(
                            MethodjoinCutThresholdSettings.
                            check_similarity_buffer).area

        ## !save 保存中间结果
        if(check_save_option):
            df.to_file(f"{self.dataprovider.check_path}/check_rp_{ldbm}.shp", encoding='utf-8')
            df_rs.to_file(f"{self.dataprovider.check_path}/check_rs_{ldbm}.shp", encoding='utf-8')
            
        
        ## 绘图
        if (image_save_option):
            fig, axs = plt.subplots(nrows=1, ncols=5, figsize=(10, 8))
            ax = axs[0]
            df_rs.plot(column="NID", ax=ax, cmap='gist_rainbow')
            df.plot(column="NID", ax=ax)
            ax.set_title("rp and rs")

            ax = axs[1]
            df_rs.plot(column="NID", ax=ax)
            ax.set_title("rs")
            # self.viz.matplot(plt)
        print("切分遥感路段", len(df_rs))
        print(f"切分上报路段{ldbm}", len(df))


        ## * 5.切割结果移交与整理  
        df["method"] = ["rp"] * len(df)
        index_rp_geom = list(df.columns).index('geom')
        index_rp_method = list(df.columns).index('method')
        index_rp_NID_rs = list(df.columns).index('NID_rs')
        for line_i in range(len(df)):
            start = df.iloc[line_i].start
            end = df.iloc[line_i].end
            x = df_rs.start == start
            y = df_rs.end == end
            logical_and = np.logical_and(x, y)
            rs_match = df_rs[logical_and]
            if (len(rs_match) != 0):
                if (bufferPer[rs_match.NID.values[0], line_i] >
                        self.THRESHOLD_BUFFER_COVER):
                    df.iloc[line_i, index_rp_geom] = rs_match.geom.values[0]
                    df.iloc[line_i, index_rp_method] = 'rs'
                    df.iloc[line_i, index_rp_NID_rs] = rs_match.NID.values[0]
        df_before_connect = copy.deepcopy(df)

        ## 获取未融合路段
        join_rs_NID  = list(df[df.method == 'rs'].NID_rs)
        df_not_join = df_rs.query("NID not in @join_rs_NID")
        
        ## * 6.连接与捕捉
        index_geom = list(df.columns).index('geom')
        while (len(df) > 1):
            geo1 = df.geom.iloc[0]
            geo2 = df.geom.iloc[1]
            if (df.method.iloc[0] == 'rp'):
                # 直接合并
                ## 因为移交了rs的，做了一些变换，所以会产生间隙，由于是rs捕捉到rp，rp本身应该没有变，rs-rs可能会有问题，rp-rp后面再看看
                if (df.method.iloc[1] == 'rp'):
                    geo_merge = shapely.ops.linemerge([geo1, geo2])
                    assert type(geo_merge) == shapely.geometry.LineString
                    df.iloc[1, index_geom] = geo_merge
                    df.drop(index=[0], inplace=True)
                    df.reset_index(drop=True, inplace=True)
                    pass
                else:  #rs
                    if (joinMethod == 'shapely'):
                        geo_snap = shapely.ops.snap(
                            geo2, geo1, tolerance=self.THRESHOLD_TORLERANCE)
                    elif joinMethod == 'fix':
                        geo_snap = self.my_extend(
                            geo2, geo1, tolerance=self.THRESHOLD_TORLERANCE)
                    elif joinMethod == 'move':
                        geo_snap = move_extend.move_snap(geo2,
                                                               geo1,
                                                               threshold=500)
                    else:
                        raise NotImplementedError
                    geo_merge = shapely.ops.linemerge([geo_snap, geo1])
                    # geo_merge = geo_snap
                    assert type(geo_merge) == shapely.geometry.LineString
                    df.iloc[1, index_geom] = geo_merge
                    df.drop(index=[0], inplace=True)
                    df.reset_index(drop=True, inplace=True)
                    pass
            else:  # rs
                if (True):  #df.method.iloc[1] == 'rp'):
                    if (joinMethod == 'shapely'):
                        geo_snap = shapely.ops.snap(
                            geo1, geo2, tolerance=self.THRESHOLD_TORLERANCE)
                    elif joinMethod == 'fix':
                        geo_snap = self.my_extend(
                            geo1, geo2, tolerance=self.THRESHOLD_TORLERANCE)
                    elif joinMethod == 'move':
                        geo_snap = move_extend.move_snap(geo1,
                                                               geo2,
                                                               threshold=500)
                    else:
                        raise NotImplementedError
                    geo_merge = shapely.ops.linemerge([geo_snap, geo2])
                    assert type(geo_merge) == shapely.geometry.LineString
                    df.iloc[1, index_geom] = geo_merge
                    df.drop(index=[0], inplace=True)
                    df.reset_index(drop=True, inplace=True)
                    pass
                else:  #rs
                    # geo_merge = shapely.ops.linemerge([geo1, geo2])
                    # assert type(geo_merge) == shapely.geometry.LineString
                    # df.iloc[1, index_geom] = geo_merge
                    # df.drop(index=[0], inplace=True)
                    # df.reset_index(drop=True, inplace=True)
                    pass
                pass
            pass  # while
        ## * 捕捉到两端
        geo_final = df.iloc[0].geom
        geo_final = move_extend.move_snap_both_end(geo_final,
                                                         geom_rp,
                                                         threshold=500)
        df.iloc[0, index_geom] = geo_final

        ## * 桩号分配与属性转移 --- 好像不需要在这里做，在外面一起了 这里废了

        if (image_save_option):
            ax = axs[2]
            df_before_connect.plot(column="method", ax=ax, cmap='spring')
            df_before_connect.to_file(f"{self.dataprovider.check_path}/roads_before_connect_{ldbm}.shp")
            ax.set_title("join")
            ax = axs[3]
            self.rp_df[self.rp_df[self.dataprovider.rp_layer.pk] == ldbm].plot(ax=ax)
            ax.set_title("origin")
            ax = axs[4]
            df.plot(ax=ax)
            ax.set_title("final")
            plt.suptitle(ldbm)
            # self.viz.matplot(plt)
            plt.savefig(f"{image_path}/{ldbm}_process.png")
        return df,df_not_join

    def joinMatchCutWorkflow_ldbmList(self,
                                      ldbmList,
                                      joinMethod='shapely'):
        """
        对满足条件2的上报路段工作流程
        备注：条件2---一个上报路段对应的遥感路段，可以连接为一整个路段
        
        Args:
            ldbmList (list): 待处理上报路网路段编码
            joinMethod (str, optional): 捕捉合并方法，默认为shapely，可选fix与move. Defaults to 'shapely'.
            save (bool, optional): 结果保存选项，默认保存结果. Defaults to True.

        Returns:
            [type]: [description]
        """
        status = []
        matchResult = []
        for ldbm in ldbmList:
            try:
                tmp_match_rodas_df = self.match_df[self.match_df[self.dataprovider.rs_layer.match_road] ==
                                                   ldbm]
                tmp_rp_rodas_df = self.rp_df[self.rp_df[self.dataprovider.rp_layer.pk] == ldbm]
                geom_rp = tmp_rp_rodas_df[self.dataprovider.rp_layer.geom].values[0]
                geom_rs = tmp_match_rodas_df[self.dataprovider.rs_layer.geom].values
                matchdf,not_match_df = self.joinMatchCutWorkflow_single(
                    ldbm, geom_rs, geom_rp, joinMethod)

                status.append("ok")
                matchResult.append({
                    "ldbm": ldbm,
                    "geom": matchdf.geom.values[0]
                })
                
            except AssertionError as assert_failure:
                print(str(assert_failure))
                status.append("assert")
            except SnapException as snap_error:
                print(str(snap_error))
                status.append("snap")
                # self.viz.matplot(plt)
                plt.savefig(f"{check_path}/{ldbm}_snap_failure.png")
            except ExtendException as extend_error:
                print(str(extend_error))
                status.append("extend")
                # self.viz.matplot(plt)
                plt.savefig(f"{check_path}/{ldbm}_extend_failure.png")
                pass
            except Exception as e:
                print(str(e))
                status.append("other")
        gdf = gpd.GeoDataFrame(matchResult,
                               geometry=self.dataprovider.rp_layer.geom,
                               crs=f'EPSG:{self.SRID}')
        if result_save_option:
            gdf.to_file(f"{self.dataprovider.result_path}/JoinResult_{self.dataprovider.rp_layer.name}.shp")
        df = pd.DataFrame(ldbmList)
        df["status"] = status
        print(df)
        return gdf,df

    def joinMatchDisconnectWorkflow_ldbmList(
            self,
            ldbmList,
            joinMethod='shapely',
            THRESHOLD_SHORTITEM_CLEAR_LENGTH = 20,
            THRESHOLD_SHORTITEM_CLEAR_COORDS_NUM = 5):
        
        """
        对满足条件3的上报路段工作流程
        备注：条件3---一个上报路段对应的遥感路段，无法连接为一整个路段，由至少两个部分组成
        
        Args:
            ldbmList (list): 待处理上报路网路段编码
            joinMethod (str, optional): 捕捉合并方法，默认为shapely，可选fix与move. Defaults to 'shapely'.
            THRESHOLD_SHORTITEM_CLEAR_LENGTH (int, optional) : 碎屑路段长度阈值，小于此长度的碎屑路段会被排除，默认为20(m)
            THRESHOLD_SHORTITEM_CLEAR_COORDS_NUM (int, optional) : 碎屑路段内部结点数量阈值，内部结点数量少于此阈值的碎屑路段会被排除，默认为5个
            layer (str, optional): 处理图层名. Defaults to "".
            save (bool, optional): 结果保存选项，默认保存结果. Defaults to True.

        Returns:
            [type]: [description]

        """
        status = []
        matchResult = []
        for ldbm in ldbmList:
            try:
                tmp_match_rodas_df = self.match_df[self.match_df[self.dataprovider.rs_layer.match_road] ==
                                                   ldbm]
                tmp_rp_rodas_df = self.rp_df[self.rp_df[self.dataprovider.rp_layer.pk] == ldbm]
                geom_rp = tmp_rp_rodas_df[self.dataprovider.rp_layer.geom] .values[0]
                geom_rs_list = list(
                    shapely.ops.linemerge(tmp_match_rodas_df[self.dataprovider.rs_layer.geom] .values))

                ## * 目前是把几个部分合起来 然后直接取出各部分
                ## * 如果说为了部件之间仍保持断开，可能还要识别一下几个部分，比如用网络做做拓扑
                ## 排除过短的路段，减少碎屑的影响
                not_match_df = None
                for geom_rs in geom_rs_list:
                    if geom_rs.length < THRESHOLD_SHORTITEM_CLEAR_LENGTH or len(
                            geom_rs.coords) <= THRESHOLD_SHORTITEM_CLEAR_COORDS_NUM:
                        continue
                    matchdf,not_match_df = self.joinMatchCutWorkflow_single(
                        ldbm, [geom_rs], geom_rp, joinMethod )  ## TODO 修改部件的问题
                    
                    if (not_match_df is not None):
                        not_match_df = not_match_df.append(not_match_df)
                    else:
                        not_match_df = not_match_df
                    geom_rp = matchdf.geom.values[0]
                status.append("ok")
                matchResult.append({
                    "ldbm": ldbm,
                    "geom": geom_rp
                })
            except AssertionError as assert_failure:
                print(str(assert_failure))
                status.append("assert")
            except SnapException as snap_error:
                print(str(snap_error))
                status.append("snap")
                # self.viz.matplot(plt)
                plt.savefig(f"{check_path}/{ldbm}_snap_failure.png")
            except ExtendException as extend_error:
                print(str(extend_error))
                status.append("extend")
                # self.viz.matplot(plt)
                plt.savefig(f"{check_path}/{ldbm}_extend_failure.png")
                pass
            except Exception as e:
                print(str(e))
                status.append("other")
        gdf = gpd.GeoDataFrame(matchResult,
                               geometry="geom",
                               crs=f'EPSG:{self.SRID}')
        if result_save_option:
            gdf.to_file(f"{self.dataprovider.result_path}/JoinResult_con3_{self.dataprovider.rp_layer.name}.shp")
        df = pd.DataFrame(ldbmList)
        df["status"] = status
        print(df)
        return gdf,df